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ABSTRACT 
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QyiJ±l 


Global  approximation  methods  for  solving  the  Abel  integral  equation 

/X  — y(s)ds  = f (x)  , 0 < a < 1 , x >_  0 , 

0 (x-s) a 

by  means  of  splines  with  full  continuity  are  considered.  The  methods  jure  based 

on  using  the  differentiated  form  of  the  above  equation.  It  is  shown  that  the 

use  of  linear  splines  in  C leads  to  a 2-a  method  for  0 < a < 1 and  the 

use  of  quadratic  splines  in  C1  leads  to  a 3-a  method,  which  computational 

experiments  indicate  is  stable  for  0 < a < 1 , though  this  is  proved  only  for 
• i,n3 

0.415  (=  2 - ) £ a < 1 • The  same  technique  applied  to  cubic  and  higher- 

order  splines  gives  rise  to  divergent  methods. 


AMS(MOS)  Subject  Classification:  45E10,  45L10,  65R05 

Key  Words:  Abel  integral  equation,  Global  approximation.  Linear  splines, 
Quadratic  splines.  Cubic  splines.  Asymptotic  error  estimates. 
Higher  order  method. 
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Significance  and  Explanation 

Abel-type  integral  equations  (see  Abstract)  occur  very  often  in  applications. 
Typical  examples  are  the  solution  of  the  falling  particle  problem,  measurement  in 
axially  symmetric  systems,  and  the  analysis  of  Brownian  motion  and  diffusion  proc- 
esses such  as  heat -conduction. 

In  this  paper,  a higher  order  global  approximation  method  has  been  developed 
for  the  solution  of  a class  of  Abel  integral  equations  (see  Abstract) , using  quad- 
ratic splines  with  continuous  first  derivative.  The  method  is  based  on  using  the 
differentiated  form  of  the  given  equation  alone.  It  is  self-starting  and  the  step- 
size  can  be  changed  at  any  knot  of  the  spline  without  added  complication.  At  each 
step  only  one  equation  has  to  be  solved,  as  compared  with  other  higher  order  methods 
which  require  solution  of  a system  of  equations  whose  coefficients  usually  involve 
a much  larger  number  of  integrals  to  be  evaluated.  The  computational  effort  re- 
quired by  our  method  is  only  marginally  greater  than  that  required  for  a linear 
spline  solution  of  the  original  equation.  Convergence  is  obtained  not  only  for  the 
approximate  solution  but  also  for  its  first  two  derivatives.  Thus  the  method  is 
economical  when  the  values  of  the  solution  and  its  derivatives  are  required  at  a 
large  number  of  points  where  usual  discrete  methods  of  computation  will  be  time 
consuming.  We  have  also  derived  very  simple  asymptotic  error  formulas  for  the 
first  and  second  derivatives  of  our  approximate  solution. 


SPLINE  APPROXIMATION  TO  THE  SOLUTION 
OF  A CLASS  OF  ABEL  INTEGRAL  EgUATIONS 


Hing-Sun  Hung 

1.  Introduction.  In  this  paper  we  consider  the  Abel  integral  equation 

(1.1)  /*  — - — y(s)ds  - f(x)  , 0 < a < 1 , 

0 Cx-s)01 

which  is  to  be  solved  for  y(x)  in  0 <x<  l.  This  is  a classical  equation,  and  it  arises 
in  many  mathematical  and  physical  problems  (see,  for  example.  Noble  [8]  and  Anderssen  [1]). 
If  it  is  assumed  that  f(xj  is  differentiable,  then  the  solution  of  (1.1)  is  explicitly 
given  by  [9] 


(1.2) 


y(x) 


sin  air 
air 


. d rx  f (s) 

1 dx  i , .1-a 

0 (x-s) 


ds) 


sin  air  . f (0)  rx 
- 1 1-a  1 


f'(s) 
.1-a 


ds]  . 


x - 0 (x-s) 

The  literature  on  the  numerical  analysis  of  Equation  (1.1)  and  related  equations  is  not 
only  extensive  but  also  contains  a great  variety  of  mathematically  independent  proposals 
for  their  solution  numerically.  A bibliography  can  be  found  in  [1] . 

Recently,  Equation  (1.1)  has  been  considered  by  Weiss  and  Anderssen  111]  and  by  Weiss 
(10],  who  applied  product  integration  techniques  to  generate  approximate  values  for  y(x) 
at  a set  of  discrete  points  by  means  of  piecewise  constant  functions  and  linear  splines. 
Benson  [2]  uses  the  same  methods  to  a more  general  equation,  proves  convergence  and  obtains 
asymptotic  error  estimates.  Their  approach,  however,  does  not  allow  the  order  of  the  method 
to  go  beyond  two  theoretically.  Brunner  [3]  examines  a direct  method  for  solving  Equation 
(1.1)  which  makes  use  of  spline  functions  where  the  usual  continuity  requirements  are  some- 
what relaxed.  To  be  precise,  he  uses  piecewise  polynomials  of  a given  degree  and  of  class 
C to  generate  approximate  solutions  for  Equation  (1.1),  and  shows  that  convergence  of  order 
m is  obtained  if  mth  degree  piecewise  polynomials  are  used.  His  method  is  basically  a 
block  by  block  method  in  the  sense  of  (7]  which  requires  at  each  step  the  solution  of  a 
system  of  equations. 

The  aim  of  the  present  paper  is  to  use  our  concepts  developed  in  [5]  to  obtain  a 

,1 


global  approximate  solution  for  Equation  (1.1)  by  quadratic  splines  in  C 


The  method 


Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-75-C-0024. 


to  be  described  in  Section  2 is  based  on  using  the  differentiated  form  of  (1.1)  alone. 

It  is  a self-starting  step  by  step  method;  the  stepsize  can  be  changed  at  any  knot  of  the 
spline  without  added  complication.  At  each  step,  only  one  equation  has  to  be  solved.  The 
computational  effort  required  is  only  marginally  greater  than  that  required  for  a linear 
spline  solution  of  the  original  Equation  (1.1)  (see,  for  instance,  [2]  or  [10]).  In 
Section  3,  we  first  derive  asymptotic  error  estimates  for  the  first  and  second  derivatives 
of  the  approximate  solution,  then  we  prove  that  the  orders  of  convergence  for  the  approxi- 
mate solution  as  well  as  its  first  two  derivatives  at  each  point  in  the  interval  of  inte- 
gration are  3-a,  2 and  1 , respectively,  for  0.415(  = 2 ~ ^ 2 > — “ < • This  interval 

contains  the  important  case  a * — . Numerical  examples  are  given  in  Section  4.  Finally, 
we  show  by  an  example,  in  Section  5,  that  the  same  technique  applied  to  cubic  splines  in 
C2  leads  to  a divergent  method.  A 2-a  method  for  solving  (1.1)  by  linear  splines  is 
also  considered.  The  results  are  summarized  in  the  appendix  of  this  paper. 

In  the  present  paper  it  is  assumed  that  0 £x  < 1 , but  this  restriction  is  not 


essential. 


. 


2.  The  Quadratic  Spline  Method.  Let  «i  ■ ih  , i ■ 0,1,2,...,  where  h la  an  arbitrary 
constant  stepsize.  Let  y(x)  be  the  exact  solution  of  (1.1),  and  let  Y^Y^  denote 
approximations  to  y(xj,  y'  (xj  , respectively.  We  use  a quadratic  spline  P(x)  in 
C1[0,1]  with  knots  at  the  points  x^  as  an  approximation  to  y(x) , i.e.,  for  i « 0,1,2,.. 

(2.1)  P(x)  “ + J t (2u^(x)  - u^(x))Yj  + u|(x)Y^+1),  x < xi+1  , 

where  uAx)  » (x-x^/h.  Differentiating  (2.1)  we  obtain  the  derivative  of  P(x)  , 

(2.2)  P’(x)  = ^(xjY^  - ui+1(x)Yj.  x.  <x  < xi+1  . 


Both  P(x)  and  P’(x)  are  continuous  at  the  knots. 

In  order  to  obtain  a numerical  solution  of  (1.1)  by  using  the  quadratic  spline  (2.1), 
we  have  to  introduce  the  differentiated  form  of  (1.1).  Assume  that  we  can  differentiate 


(1.1) . Then 


0 <x-s)1 


y'lslds  - f'(x)  - y(0)x  . 


The  approximate  solution  P(x)  of  the  integral  equation  (1.1)  is  derived  from  values 


P(x.)  = Y.  and  P'(x.)  = Y!  . We  know  how  to  find  Y'  ,,  assuming  that  we  know  Y!  , 
li  li  k+1  ’ i 

i = 0,1,..., k.  We  can  then  deduce  Y^+^  from  the  equation  obtained  by  setting  x = 
in  (2.1),  which  gives 

(2.4)  Y = Y + — (V  + Y'  ) 

11  k+1  k 2 v'k  k+V  - 

To  compute  Y^  , we  require  that  P’ (x)  satisfies  (2.3),  i.e.,  y'(x)  is  replaced  by 

P*  (x)  derived  from  the  values  P'lxA  = Y|  computed  from 
x 

(2.5)  /k  — ^ P'  (s)ds  = f ' (x  ) - y(0)x~a  , k = 1,2 

0 (xk-s)a 


This  can  be  rewritten  in  the  form: 
k 

(2.6)  T w,  ,y:  = f'(xj 


lo  wk.iYi  ‘ ' 


k - 1,2 


1 I 


3.  Convergence  Results.  The  proof  of  Theorem  3.1  requires  the  following  lemna  which  is 
a consequence  of  the  standard  results  for  regular  infinite  system  of  algebraic  equations 
by  Kantorovich  and  Krylov  in  16,  p.  27] . 


I 


Lemma  3.1.  If  there  exist  a constant  c > 0 and  an  integer  N > 1,  all  independent  of 


k , such  that 


|xj  ±c  . 


Vll  i l K+l,iHXil  + iBj' 


with 


i-0 
k 


i - 0,1,. ...N 

k - N ,N+1 , . . . , 


pk  * 1 ' J0  lvi.il  > 0 ' 


le.  I < c p , 


k = N ,N+1 , . , 


then 


l*£l  i c » 


i - 0,1,2,. 


The  proof  of  this  lemna  can  be  found  in  [4,  p.  3] . 

Let  y(x)  be  the  exact  solution  of  (1.1)  and  define  the  discretization  error  function 


e(x)  by  e(x)  = y(x)  - P(x),  where  P(x)  is  the  quadratic  spline  approximation  to  y(x) 


obtained  frcm  our  numerical  method.  Denote  (xj  by  c!r^  (r  = 0,1,2).  We  state 


the  following  theorem  for  the  asymptotic  estimate  of  : 


Theorem  3.1.  If  y e C4[0,1]  and  2 - <_  a < 1 , then  for  k = 1,2,. . . 


(3.1) 


ek  ° II  y”  (xk}  + 0<h2A1_a> 


4 

Proof:  Since  by  assumption  y(x)  € c [0,1],  it  is  not  difficult  to  show  that,  for 


x € Ixi'.*i+1]  , 


(3.2) 


e(x)  = ei  + j [(2ui(x)  - u3  (x) ) e|  + u3(x)  e^+1) 


h 3 5 

+ -JJ  y'"  (xA)  (2ui(x)  - 3u^(x))  + p(x)  , 


where 

(3.3) 


p(x)  = n 12  i*  > V4)  (s)ds  - 3u2i(x)  /*i+1  vi,8l»141  (s,dsI  ' 


with  u^(x)  as  defined  in  (2.1), 


x . 
l 


Adding  and  subtracting 


(x)y,M  (x.)]  to  the  right  side  of 


(3.4)  and  using  the  Taylor  series  expansion  for  y'" (x)  at  x 


difficult  to  show  that 


Since  both 


(i+1) 


Now  multiply  (3.10)  by  k 

(k+i)V  , 


divide  by 


to  yield  the  required  error  equation 


(k+2) 


with 


y'" (x  ) = 0(h  ) , it  is  easily  shown  from  Equation  (3.10)  by  using 


with  K as  defined  in  Lemma  3.2 


below.  On  the  basis  of  this  together  with  Lemma  3.2(b)  and  Lemma  3.3(c)  below  we  can  apply 


Thus  the  proof  of  Theorem  3.1  is  completed 


Theorem  3.2.  Let  the  assumptions  of  Theorem  3.1  be  satisfied,  then  for  lc  * 1,2 


in  (3.2)  and  noting  that  p(x.  ) * 0(h  ),  we  obtain 


0 , by  means  of  (3.1)  and  using  the  fact  that  J, 

i=l 

(3.17)  follows  immediately  from  (3.18)  by  induction 


ym  (x  ) + 0(h/(k+l)  A-a) 


Proof : Differentiating  (3.2)  twice  and  letting  i = k we  have  for  x < x < x 


Setting  x = x in  (3.20)  and  using  (3.1),  (3.19)  immediately  follows 


■■  - --  - ■.  ^ - V-  . 


In  the  following  we  can  show  that  the  approximate  solution  P(x)  along  with  its 


first  and  second  derivatives  converges  to  the  corresponding  exact  solutions  at  each  point  in 
the  interval  of  integration. 

Theorem  3.4.  Let  the  assumptions  of  Theorem  3.1  be  satisfied,  then  for  any  fixed  x e [0,1] 
| e(x)  | = 0(h3  “)  , | e'  (x)  | = 0(h2)  , | e"  (x)  | = 0(h)  . 

Proof : By  means  of  (3.1)  and  (3.17),  Theorem  3.4  follows  immediately  from  Equations  (3.2), 

(3.4)  and  (3.20). 


The  result  of  Theorem  3.4  has  been  obtained  under  the  assumptions  that  y e C [0,1] 

and  that  2 - ^ ^ <.  <*  < 1 • If  we  relax  the  assumption  on  y , viz,  y'"  is  Lipschitz 

continuous  in  [0,1] , we  can  show  that  convergence  for  the  approximate  solution  and  its  first 

two  derivatives  is  still  possible,  the  order  of  convergence  being  2,  2 and  1 , respectively. 

We  can  also  relax  the  assumption  on  a by  modifying  the  structure  of  the  a,  , . 's  in  the 

k+l,i 

error  equation  (3.12)  using  a technique  analogous  to  that  of  [10].  But  it  will  involve  a 
complicated  analysis.  Since  our  proof  already  includes  the  important  case  a = j , an 
extension  of  Theorem  3.4  in  this  direction  will  not  be  pursued  here.  Numerical  experiment  in 
Section  4 indicates  that  the  result  of  Theorem  3.4  holds  for  all  a , o < a < 1 . 


In  3 


< a < 1 , then  there 


Lemma  3.2.  If  the  a . 's  are  defined  by  (3.13)  and  2 , 

K+  l/i  Jen  2 — 

exist  integers  K and  K >_  1 , independent  of  h and  k , dependent  on  a , such  that 


(a) 


Vi.i  - 0 ' 

•kti.ki0 ' 


i = 0,1,...  ,k-l,  k >_  1 , 
k > K , 


(b) 


IJVl.il  > (2-a,(1-a> 


i=0 


(1— c) 
(k+l)° 


for  an  appropriate  c , 


a < c < 1 


k > K , 


Proof  of  (a) : From  (3.15),  using  (2.7),  we  have 

x. 


I1 


k+1'°  (k+l)aw  x 

k+l,k+l  0 


*1  (Xj-s) 


(1-lf 

■ Xk 


U-x— >“ 

Vi  j 


l+a 


ds  > 0 


since  the  integrand  is  nonnegative.  Similarly,  we  can  prove  that  i — ® f °r  i = 

1,  ...»  k-1  . Finally,  it  is  easy  to  show  that  as  k increases  a,  , , tends  to  3 - 

k+1  ,k 

which  is  greater  than  or  equal  to  zero  for  2 - "j— 3 £ a < 1 , therefore  there  exists  a 

integer  K , independent  of  h and  k such  that  a,  , v > 0 for  k > K and 

k+ 1 » K — — 


,2-“ 


2 - . - -j  < a < 1 . Thus,  part  (a)  follows  immediately  from  (3.13). 

k 

Proof  of  (b) : Since  it  would  be  complex  if  we  estimate  J a . directly,  we  introduce 

i=0  (1 


k+l,i  ‘ 


(3.21) 


,k+c  . l-o 

k k,i  k+l,i 

Wk+l,k+l 


, 0 < c < 1 , i=0,l,...,k 


- . ,k+c.l-a  , k .a  _ , , ...  . . , tn  3 

By  noting  that  (— — j > (: — — ) for  k > 1 , it  is  obvious  that  for  2 - - — — < a < 1 , 

k k+1  — In  2 - 

(3-22)  = °'1 k'1'ki1' 

•kfi.k  ” Vi.k  i 0 - * i K ' 

where  the  K is  the  same  K as  defined  in  part  (a)  . 

Since  from  (2.7)  we  have 

k x.  xv  a 

(3-23)  l wk  i = / — —-‘fc-irs- 

i=0  R,i  0 (x.-s)a  1 “ 

k 

for  k = 1,2,...  , it  can  easily  be  shown  that  for  2 - — ^ a < 1 , 


1 - l l«k+lfil  1 (2 -a)  (1-a) 


k > K . 


If  we  can  show  that  for  an  appropriate  c , 0 <_  c < 1 , a - a . . 0 , then  the 

proof  is  completed.  To  do  this,  it  is  sufficient  to  show  that,  for  i = 0,1, ...,k 

_ k+c  .l-o,  k+1  .a  , , ,,  a . k+2  .1-0  , ,,,  a ,,  ,,a  . 

Dk+i,i=[(  — > ( — » - 11[kwk,i]  ' t(  - mkWki  - (k+D  w^ij  >° 

By  using  (2.7)  and  letting  s = ht,  it  is  easily  verified  that  for  i = l,...,k-3, 

(3.25)  Vl,i  - t/i  Lk  t (t-i.+l)dt  + /1+1  Lk  i(t)  (i+l-t)dt]h1_“  , 


L.  . (t)  = [(  ^■)1'°‘(  ^±1)“  - 1)  — [(  ^f)1_a  - 1]  -5—^ . 

k'i  k k (k-t) a 1 k1  “(k-t)1'01 

Since  for  i-1  _<  t _<  i+1,  i = l,...,k-3  and  o <_  c < 1 , 

, 2 

L.  . (t)  > 1-2-  > 0 
k,i  - k2 

when  k is  sufficiently  large,  it  is  obvious  from  (3.25)  that  i — ^ for  i “ 

2.n  3 

1,  ...,  k-3  . Also  when  a c < 1 and  2 - ^ ^ <_  a < 1 , it  is  easy  to  show  by  using 


y ft 


i 

t; 


i : 


Let  M - max  |y'"  (x)  | and  M « max  |y  4 (x)  | . Then  by  straightforward 

J r a 11  ' m 


xt [0.1] 


xf [0,1] 


estimation  and  noting  that  hk  <_  1 , we  obtain  from  (3.27) 
(3.28) 


£ | , allta)  , 1 1 1 . , v3-'1  . , 5 ,4-a  1-a 

Rk'  - 1 24  1-a  16  1+a*1  M3h  1 12  (1-a) 5 M4h  k 


— ci  h' 


3-a 


where  = [ 


a (1+a) 


. 1 11,, 

(tt  + — rrr  > 1 m„  + 


24  '1-a  16  1+a  ' J "3  T 12  (1-a)  4 * 

Proof  of  (b)  : Subtraction  of  (3.27)  from  (3.27)  with  k replaced  by  k+1,  then  by 

straightforward  estimation  and  noting  that  hk  £ 1 , it  is  not  difficult  to  show  that 
for  k = 1,2, .. . 

.-(1)  .*(1),  . i - ( 2)  - C2)  | 


(3.29) 


IVi-Nlil^i-*k  l+  I^i-^l 


. ,o(l+o),  1 .1  1 .4-a  . a (1+a)  1 ..  h3"“  , 

< t ^ iTir  + rr  rrri  M^h  ♦ — 72 — — **,  -rrr  ) 


- ' 24  '1-a  16  1+a'  4 


+ (1  + ifc  > M4h4"“  > 


24  16  3 ^2+a 


3-a 


- h 
± °2  — 


where  c„ 


ad+a)  M . / 

384  M3  { 24 


Proof  Of  (c) ; From  (3.14),  using  (3.28),  (3.29)  and  (2.7),  we  obtain  for  k = 1,2,. 
|b  | < ,(lc+2)1"--  - I ■ ” - k 


k+1, k+1 

< *3  ^ 

“ 3 k“ 


{lRk+i  - + 11  “ ( k7T>aH\l> 


where  c = 3 (1-a)  (2-a)  (a  c +c  ) 
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5.  Divergence  of  Higher  Order  Spline  Methods.  We  have  shown  that  the  linear  spline 

method  (see  Appendix)  and  the  quadratic  spline  method  (.see  Section  2)  based  on  using  the 

differentiated  form  of  Equation  (1.1)  alone  are  convergent.  It  might  be  supposed  that  the 

same  technique  applied  to  cubic  and  higher  order  splines  would  result  in  higher  order 

approximation  methods  for  Equation  (1.1).  Such  an  approach,  however,  fails,  because  of 

instability,  in  the  case  that  the  spline  is  of  full  continuity.  For  illustration,  apply  the 

2 

method  by  using  a cubic  spline  in  C [0,1]  to  the  following  Abel  integral  equation: 

rx  1 6 4- a 

0 (x-s)  ° y(S  dS  * U-a)  (2-a)  (3-a)  (4-a)  X ' 

where  y(x)  =■  x3,  with  a = 0.1,  0.5  and  0.9.  The  divergence  of  the  numerical  results 

can  be  seen  very  clearly  in  Table  5.1. 
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APPENDIX 


Solving  a Class  of.  Abel  Integral  Equations 
By  Linear  Splines 


Introduction.  In  this  appendix  we  consider  the  use  of  a linear  spline  to  obtain  a 
global  approximate  solution  for  Equation  (1.1).  We  shall  describe  the  method  briefly 
and  state  the  important  theorems.  The  proofs  of  the  theorems  have  been  carried  through 
but  the  details  are  emitted,  since  the  methods  of  proof  are  similar  to  those  already  used 
in  this  paper. 

A Linear  Spline  Method.  Let  x^  = ih,  i > 0,1,...,  where  h is  an  arbitrary  constant 

stepsize.  Let  Y^,Y|  denote  approximations  to  y(x  J , y ' (xj  , respectively.  We  use  a 
linear  spline  function  P(x)  in  CI0,1]  with  knots  at  the  points  x^  as  an  approxima- 
tion to  y(x),  i.e.,  for  i = 0,1,2,... 

(A.l)  P(x)  = Y,  + Y!(x-x.),  x.  < x < X.,,  . 

i l i i — — x+1 

The  function  P(x)  is  continuous  at  the  knots. 

The  approximate  solution  p(x)  of  the  integral  Equation  (1.1)  is  derived  from  values 


P(xJ 


Yi  and  P'ixJ 


y;  . We  know  how  to  find  Y'  , assuming  that  we  know  Y'  , 

X K X 


0,1,..., k-l.  We  cam  then  deduce  Y 


k+1 


from  the  equation  obtained  by  setting  x = x 


in  (A.l),  which  gives 
(A.  2)  Y, 


‘k+1 


Yk  + h Yk  • 


To  compute  Y|  , we  require  that  P'(x)  satisfy  (2.3),  i.e. 
derived  from  the  values  P ' (x^)  =■  Y|  computed  from 


k+1 


y'(x)  is  replaced  by  P'(x) 


(A. 3) 


k+1 


lxk+rs) 


- P'  (s)ds  = f • (x^?  - y(0)xfc+1  , k = 0,1,2,... 


This  can  be  rewritten  in  the  form: 
k x. 


(A. 4) 


l </ 

i-0  x. 


‘i+1 


(*k+rs) 


- ds)Y-  = £'(xk+1)  - y(0)xk+1  , 


0,1,2,... 


Equation  (A. 4)  is  a nonsingular  triangular  system  for  Y^ 


For  the  starting  value  Y , we  take 


(A. 5) 


y(0)  = lim  (1-a)  — : 

1-a 

x+  0 x 
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